doseresNMA antidep/drnma.model.class.R

In htx-r/doseresNMA:

drnma.class <- function(){

  for(i in 1:ns){ # loop through ns trials
    w[i,1] <-0
    delta[i,1] <- 0
    u[i] ~ dnorm(0,0.001) # baseline effect
    for (k in 1:na[i]){ # loop through na arms within a study
      r[i,k] ~ dbin(p[i,k], n[i,k])
      logit(p[i,k]) <- theta[i,k]
      theta[i,k] <- u[i] + delta[i,k]

      # arm deviance
      rhat[i,k] <- p[i,k] * n[i,k] # the expected number of events
      resdev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (n[i,k] - r[i,k]) * (log(n[i,k] - r[i,k]) - log(n[i,k] - rhat[i,k]))) # the resuidual deviance
    }
    # study devinace
    resstudydev[i] <- sum(resdev[i, 1:na[i]])

    for(k in 2:na[i]){ # treatment effects: delta's
      # RE delta's
      delta[i,k]~dnorm(md[i,k], prec[i,k])
      # dose-effect model
      DE[i,k] <- ((beta1[class[i,k]] * dose1[i,k]) + (beta2[class[i,k]] * dose2[i,k])) - ((beta1[class[i,1]] * dose1[i,1]) + (beta2[class[i,1]] * dose2[i,1]))
      # multi-arm correction
      md[i,k] <- DE[i,k] + sw[i,k] # modify the mean
      prec[i,k] <- prec.d *2*(k-1)/k # modify the precision
      w[i,k] <- delta[i,k] - DE[i,k]
      sw[i,k] <-sum(w[i,1:(k-1)])/(k-1)
      }
  }

#** PRIORS
  # 1. to  delta's RE model
  tau.d~dunif(0,5)      # common standard deviation is given a vague prio
  prec.d<-pow(tau.d,-2) # common precision in RE-NMA model

  # 2. to beta's (DR coeff)
  beta1[refclass.index] <- 0
  beta2[refclass.index] <- 0
  for (c in c(1:(refclass.index-1),(refclass.index+1):nc)){
    beta1[c]~dnorm(0,0.001)
    beta2[c]~dnorm(0,0.001)
  }
  tau.beta~dunif(0,5)      # vague prior to beta's SE
  prec.beta<-pow(tau.beta,-2) # precision of beta's

  # overall deviance
  totresdev <- sum(resstudydev[])

  #** drug absoulte effect
   # the baseline effect
   for (m in 1:ns){
    rr[m,1] ~ dbinom(p0[m],nn[m,1])
    logit(p0[m]) <- OO[m]
    OO[m] ~ dnorm(Z, prec.z)
   }
  # baseline priors
  Z ~ dnorm(0, 0.001)
  prec.z <- pow(sigma.z,-2)
  sigma.z ~ dunif(0,10)
  p.placebo <- exp(Z)/(1+exp(Z))
  # combine baseline with OR -> drug absolute effect
   for (c in c(1:(refclass.index-1),(refclass.index+1):nc)){
    for( j in 1:nd.new){
      OR[j,c] <- exp(beta1[c]*new.dose[c,j]+ beta2[c]*f.new.dose[c,j])
      odds.drug[j,c] <- OR[j,c]*exp(Z)
      p.drug[j,c] <- odds.drug[j,c]/(1+odds.drug[j,c])
    }
  }
}